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1. Introduction 

1.1. Tensor network states 

While the density matrix renormahzation group (DMRG) [T] has proven to be a very 
useful numerical tool for the study of strongly correlated quantum many-body systems 
in one dimension, applications to higher dimensions have remained very limited so far 
because an exponential growth of the complexity with system size is recovered in general. 
Monte Carlo simulations, on the other hand, are not limited by the dimensionality of 
the problem, but perform well only for non-frustrated bosonic systems. Extensions of 
the DMRG method to two dimensions are therefore highly desired. 

In recent years, insights from quantum information theory have led to an 
understanding of the success of DMRG in terms of entanglement of the ground state: 
the renormahzation procedure of choosing the dominant eigenvectors of the reduced 
density matrix picks out states with a large contribution to the entanglement entropy. 
It can therefore be regarded as a low-entanglement ansatz for the ground state; the 
entanglement that can be encoded is bounded by aS < logM, where M denotes the 
number of states kept in a renormahzation step. Since for non-critical systems the 
amount of entanglement saturates as a function of system size [2j, a finite M will be 
able to capture the physics of such systems even in the thermodynamic limit. For weakly 
entangled systems, the ansatz will be a very accurate description even for very small M. 

Since the block entanglement does not saturate in the thermodynamic limit for 
two-dimensional systems, the matrix size needs to grow exponentially with the system 
size in an approach based on DMRG to obtain a fixed accuracy [3j. However, one can 
hope that the entanglement grows less than the volume of the system; indeed, it has 
been shown that for many classes of interesting strongly correlated systems an area law 
holds [4J, i.e. that the entanglement between a given block within a system and the rest 
of the system will scale with the surface between the two. Therefore, algorithms that 
obey an area law by construction appear as a viable approach to higher-dimensional 
systems. One example for such an algorithm will be the topic of this paper. 

DMRG can be understood as a variational method operating on matrix product 
states [5j. For states in a Hilbert space spanned by an orthonormal basis |n) = 
rii^i l(^I)^'I^)' where the a] are appropriate creation operators, n = ni . . . is an 
occupation number vector and \Q) is the vacuum, the coefficients c(n) of the wave 
function |^^) = c(n)|n) are expanded as product of a set of matrices of size M x M, 
where one matrix is associated with each lattice site: 

c{n) = Al[m]Al^^[n,]...A^Jn^] (1) 

In the limit of infinite matrix size, the set of matrix-product states is complete. It has 
been shown that MPS offer an efficient description of local and non-local properties of 
ground states of gapped and even critical systems j6]. 

For higher-dimensional systems, several ansatz states and algorithms have been 
proposed that capture an area law by construction. We will focus on an ansatz that 
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Figure 1. (color online) From matrix-product states to projected entangled-pair states. 
In both linear operator mapping from auxiliary spaces on the bonds is 

associated with each site. 

was first suggested by G. Sierra and M. A. Martin-Delgado under the name of vertex 
matrix product ansatz [7j and subsequently applied to evaluate the partition function of 
a 3D classical model by Nishino and Okunishi [El A first application to infinite, 
homogeneous quantum systems is found in |T0]. A formalism for inhomogeneous, 
finite systems was proposed by Verstraete and Cirac [llj under the name of projected 
entangled-pair states. In [12j, infinite systems were reconsidered from the point of view 
of extending Verstraete and Cirac's formalism. For this paper, we choose to refer to the 
ansatz as PEP states. 

Since it was first proposed, many studies of classical and quantum systems using 
this family of states have appeared [l3l[H[l5l[l6l[I7l[l8l[ia[20l[^ Additionally, 
other algorithms operating on these states have been proposed, including tensor 
entanglement renormalization group (TERG) E5j. Related proposals include 

the multiscale entanglement renormalization ansatz (MERA) [26j, string-bond states 
[27] and entangled-plaquette states [28] . 

1.2. Projected entangled-pair states 

In the following, we will describe PEP states, which can be understood easily from a 
diagrammatic representation (Fig. [T]). For spins distributed on an arbitrary graph, a 
tensor is associated with each vertex; the tensor has one index for each edge emanating 
from it plus one physical index (a specific choice of basis is implied). Specializing 
to the square lattice in two dimensions, we have rank 5 tensors. In the following, 
we will choose the dimensions of all indices between tensors equal to M and denote 
the dimension of the physical Hilbert space by d. The bond dimension M provides a 
systematic refinement parameter: while for ground state properties, an M = 1 PEPS 
corresponds to a mean-field state, the class of states is complete in the limit of M ^ oc. 
The values of M attainable in practical calculations are very limited due to a strong, 
but polynomial increase of the computation time with M . On the other hand, the fact 
that for increasing dimension of the lattice, typical ground states of local Hamiltonians 
are closer to mean-field states suggests that the dimension of the bonds in a PEPS can 



Assessing the accuracy of projected entangled-pair states on infinite lattices 



4 



be chosen smaller than in an MPS. Nevertheless, this is the most serious limitation to 
the accuracy of the algorithm. One of the main purposes of this paper will be to assess 
the accuracy obtained with PEPS for the values of M attained in practical calculations. 

For a wide variety of systems, the PEPS may provide a very efficient representation 
of the exponential number of coefficients with a small number of parameters. However, 
the aforementioned results on the accuracy of MPS for gapped and even critical systems 
cannot be extended to PEPS. Additionally, it can be shown that the exact calculation 
of an expectation value based on such a network, i.e. the trace over all indices, is an 
exponentially hard problem and therefore intractable directly for larger system sizes 
[29j. To overcome this limitation, approximate contraction schemes have been devised. 
We will discuss the infinite projected entangled-pair states algorithm (iPEPS), which 
operates directly in the thermodynamic limit. 

The treatment of the infinite lattice depends on the observation that the infinite 
tensor network can be regarded as transfer operator acting on some boundary state and 
all correlators can be expressed in terms of the dominant eigenvectors of this transfer 
operator. This represents a renormalization scheme that replaces expectation values 
of local operators with expectation values of effective operators that take into account 
long-range correlations, but can be evaluated locally. 

It should be emphasized that we do not handle boundary conditions explicitly. In 
some cases, the boundary state that we use to locally approximate expectation values 
becomes degenerate, in which case a choice of boundary condition is implied by the 
choice of boundary state from the manifold of degenerate states. For systems with 
a broken local symmetry, this degeneracy is related to the degeneracy of states with 
different symmetry breaking direction. We however predetermine this direction by the 
choice of the initial state, which is not symmetric and therefore favours a symmetry 
sector. In the case of topological order, this situation is much more intricate. For 
the purpose of this paper, we therefore restrict ourselves to systems where only local 
symmetries are relevant. Systems with topological order have been treated elsewhere, 

e.g. [HOIEIIESJ. 

The outline of the paper is as follows: in section [2| we will illustrate the details of 
the infinite projected entangled-pair algorithm and discuss properties of the ansatz in 
the vicinity of a first-order phase transition. In section [3| we apply the algorithm to well- 
known non-frustrated systems and assess the accuracy of the results in comparison to 
Quantum Monte Carlo simulations. In section |4| we extend the simulations to frustrated 
systems where Monte Carlo methods fail. 

2. Infinite PEPS algorithm 

2.1. Approximation scheme 

Several schemes are available for the approximation of the infinite tensor network by a 
small tensor network in which operators can be evaluated efficiently. This includes a 
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Figure 2. (color online) MPS-based scheme for the reduction of the infinite tensor 
network to a small network. In this scheme, diagonals or rows/columns of the lattice are 
regarded as infinite transfer operators and operators are evaluated using the dominant 
eigenvectors of these operators. We use a diagonal contraction scheme. In (a) and (b), 
the physical indices have already been traced out. (a) Tensor network with two-tensor 
unit cell, diagonal transfer operators and boundary IMPS, (b) The two-dimensional 
lattice is reduced to a quasi-one-dimensional system by approximation with the 
dominant eigenvectors of the transfer operators, (c) Using the dominant eigenvectors 
of the one-dimensional strip, a finite network is obtained in which operators can be 
applied to the physical indices. The trace over this network equals the expectation 
value of that operator. 

scheme based on boundary infinite matrix-product states as eigenvectors of diagonals or 
rows/columns of the lattice (cf. Fig. [2]) [ElEQlIin] and the corner-transfer matrix scheme 
originally proposed in the context of tensor-product states by Nishino and Okunishi \32\ 
and recently revisited by Orus and Vidal [2T]. 

While the direct comparison of the two methods in [2j by example of the two- 
dimensional quantum Ising model in transverse field yields slightly better results for the 
corner-transfer matrix, this is expected to be due to the fact that a very small bond 
dimension M = 2, 3 is sufficient to capture the physics of the Ising model. In this 
paper, we will consider models with more entangled ground states where errors result 
from insufficient M instead of an inaccurate contraction of the environment. Therefore, 
we do not expect significant improvement from the use of the corner-transfer matrix. 
We have checked this for M = 2 by example of the Heisenberg model (Section |3.2D 
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Figure 3. (color online) Graphical representation of eqn. ([2]); this network is identical 
to Fig. |2](c). The effective renormalized operator consists of the two-site operator and 
the environment representing the infinite lattice. By means of this renormalization, it 
takes into account long-range correlations in the state. The trace over all bonds yields 
the expectation value of O. 



and the dimerized Heisenberg model (Section 3.3) and found the deviations between 
the boundary MPS method and the CTM to be neghgible. In the case of the frustrated 
model (Section [4]), we use a slightly different scheme to update the boundary MPS in 
order to accomodate for the enlarged unit cell. The deviations in this case are slightly 
larger, but remain within about 2 % for the order parameter. 

Finding the dominant eigenvectors of the infinite transfer operator can be performed 
in several ways. Most commonly, the iTEBD algorithm for non-unitary evolution [33] is 
applied. We find that for lattices with a larger unit cell, better numerical stability can 
be achieved using an algorithm based on the optimization of fidelity without requiring 
orthogonalization; such an approach was suggested in [34j. For the boundary of the 
strip (cf. Fig. [2] (b)), it is sufficient to use a single tensor as boundary. 

In order to reduce the complexity of operations, transfer matrices are not stored 
explicitly. Instead, iterative schemes are applied where the matrix- vector product is 
implemented by exploiting the tensor structure in such a way that only two tensors are 
contracted at a time, storing the result in a temporary tensor. This operation can be 
mapped to a matrix multiplication, for which efficient computer libraries exist. The 
optimal order of contractions has to determined depending on the dimension of each 
bond. 

The local environment provides a means of expressing diagonal or off-diagonal 
operators O in an effective form acting on PEPS tensors, 

{^'\Om = iA\B'yO,s{AB), (2) 

where the (A, B) are PEPS tensors that can be regarded as vectors in C^^^^. The 
accuracy to which expectation values of effective operators approximate the exact 
expectation value is controlled by the bond dimension of the boundary matrix- 
product states, M^ The computational cost of increasing is very large, since the 
iTEBD algorithm requires the singular-value decomposition of a matrix of dimension 
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M'^M' X M'^M' , which scales as M^M^^; for reahstic choices of M' and M, where 
M' > M^, this becomes the most time-consuming operation of the algorithm. In 
the simulations, we choose large enough such that the results do not change with 
increased M'; nevertheless, we should mention that errors made at this step enter both 
into the calculation of the ground state and the evaluation of expectation values. 

In order to project the state into the ground state, imaginary time evolution 
is performed. For a Hamiltonian that can be written as a sum of strictly local 
terms, H = hij^ the operator exp{—f3H) can be expanded straightforwardly by 

discretizing time (3 = N6t and applying a Trotter decomposition [35j, 

U = exp{-5rH) = exp{-5rh,j) + 0{5r^). (3) 

{id) 

In the case of imaginary time evolution. Trotter errors do not accumulate and can largely 
be eliminated by reducing the time step during the simulation and using higher-order 
expansions. 

Applying the evolution gate directly to the PEPS tensors and transforming back to 
the previous form increases the bond dimension of the PEPS. Therefore, time evolution 
can only be approximated by finding the state |^^^) that minimizes 



D =11 U\^) - \^') II . (4) 

In order to achieve an optimal approximation, the norm || • || has to be chosen in such 
a way that the correlations in the system are taken into account. To this end, the 
evolution operator U can be renormalized to Ue^ using the scheme discussed above. 
Equation Q can then be written as 

D = {^\U^U\^) - {^'\U\^) (5) 
_(^|[/t|^/) + (^/|^) (6) 

= {ABy{U^UU{AB) (7) 

-{A,B'yU,^{A,B) (8) 

-{AByul^{A\B^) (9) 

+ {A\By\^{A,B^) (10) 



Minimizing this in A^ and B^ leads to a set of coupled linear equations which is 
solved iteratively by fixing one tensor and optimizing in the other, sweeping back and 
forth until convergence is reached. 

In the following, we have used two different unit cells: a staggered AB unit cell, 
as shown in Fig. [2| which allows for ferromagnetic and antiferromagnetic order, and 
an enlarged four-site unit cell, which allows more general types of order (cf. also Fig. 
[s]). We start each simulation with a large timestep of the order of 5r ^ 10~^, which is 
decreased to St < 10~^ during the simulation. Convergence is checked not only for the 
energy, but also for all other local observables that are being evaluated. 

As initial state for calculations with M = 2, we use either a completely random 
state or a product state, e.g. the Neel state, with a small amount of disorder added. 
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For calculations with larger M, we start from extrapolations of previous results to save 
computation time. 

2.2. iPEPS at a first- order phase transition 

The iPEPS method is very suitable to first-order phase transitions. In a previous work 
by Orus et al a phase transition in the anisotropic quantum orbital compass model 
was studied and the transition point located using adiabatic time evolution. We will 
take a diflFerent approach based on the observation that the dynamics close a first- 
order phase transition diflPer fundamentally between simulations on finite lattices and 
on infinite lattices with translational invariance. 

If the system is taken across the phase transition from phase A to phase B, 
phase A becomes meta-stable. Since phase A cannot locally be continuously deformed 
into phase B (assuming pure states), the dominant process by which the meta-stable 
state decays is the condensation of finite clusters of phase B (see, e.g.. Binder [36j), 
thereby breaking translational symmetry. This process is suppressed by enforcing 
translational symmetry. The imaginary time evolution we simulate therefore has a 
strongly suppressed probability of tunneling from one phase to the other in the phase 
coexistence regime where the energy densities of phase A and B are similar. 

We can therefore expect to find only homogeneous systems of one phase in the 
vicinity of the transition. By preparing a state deep within one phase and using 
this state as initial state for the imaginary time evolution of the Hamiltonian with 
different parameters, we can find the energy density of each phase for a relatively wide 
space of parameters. The transition point can then be located to high accuracy as the 
crossing point of the ground-state energy densities of the two phases; simultaneously, 
the existence of this stability allows us to clearly distinguish the nature of the phase 
transition. 

3. Non-frustrated test cases 

We study several non- frustrated spin-| models on the square lattice, where results can 
easily be compared to previous results and Quantum Monte Carlo calculations, which 
we perform using the ALPS package [371 1^- As a first test case and complementary 
to previous results [19], we study a first-order phase transition in the XXZ model in 
external field, which maps to the Bose-Hubbard model in the hard-core limit. We 
assess the accuracy in the case of gapless, long-range ordered states by example of the 
XY model and the isotropic Heisenberg model, which we extend to a dimerized model 
exhibiting a quantum critical point. 
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Figure 4. Schematic phase diagram of the XXZ model for J^y = 1, > 0, > 0. 
Indicated are the symmetries of the Hamiltonian and how they are broken in the 
respective phase. The phase diagram is symmetric with respect to the hne h = 
^^Jz ^ 0- The transition from the antiferromagnet to the canted phase is of first 
order and simultaneously changes the Z2 symmetry breaking from a staggered to a 
uniform pattern and breaks the U(l) symmetry. The transition from the canted phase 
to the paramagnet is of second order and restores U{1) symmetry. The dotted line 
from h = Jz = to h = 0^ Jz = I does not indicate a phase transition and is meant as 
guide to the eye. 



3.1. XXZ model in external field 

The Hamiltonian of the XXZ model for spin-| degrees of freedom is given by 

^ = E (^(^>7 + + J^^t^]) +hY.<yt (11) 

where ^a^ ^ denote Pauli matrices and a+ = a^ + ia^, a~ = — ia^ . By identifying 
= aj = = ^^2^, we find (up to a constant) the Bose-Hubbard model in the 

limit of hard-core bosons, given by 

H = -t ^(aja^ + a]a^) + V ^ UiUj + ji ^ n^, (12) 

where we use Jxy = —^^ Jz = j^h = . 

The phase diagram of this model has been studied to much detail and is well 
understood [42l HHJ |39]. A sketch along with the symmetries of the Hamiltonian and 
the ground state in each phase is shown in Fig. [4j In the following discussion, XY 
and Z identify the axes in spin space, not in real space. Three different phases can be 
identified: 

Paramagnet A large magnetic field forces the system into a paramagnetic state, 
corresponding to a full or empty state in the bosonic picture. The U{1) symmetry 
in the XY plane is not broken. 
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Figure 5. Spin-flop transition of the XXZ model. Lines are guides to the eye and 
interpolate the results for an iPEPS with bond dimension 2. The QMC results are 
taken from Yunoki ^39j . 
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Figure 6. Convergence of relative errors Ax = |x(iPEPS) - x(QMC)|/|x(QMC)| with 
bond dimension M. Monte Carlo results from [4^ |41J are taken as exact. For the 
deflnition of the magnetization refer to the text. 

Canted phase This state, referred to as superfluid in the bosonic language, 
spontaneously breaks the U{1) symmetry. At the line /i = 0, the Z2 Ising symmetry 
is restored. At the isotropic point Jxy = J^, this symmetry is enlarged to a U{1) 
symmetry. 

Antiferromagnet The antiferromagnetic phase (checkerboard solid) is characterized 
by a different symmetry breaking pattern of the Z2 symmetry: instead of the order 
induced by the external field, an antiferromagnetic pattern emerges. The U{1) 
symmetry is not broken. 

Only the superfluid phase is gapless. 

The transition from the full or empty solid to the superfluid phase is a second order 
phase transition which can be located analytically as the appearance of a single boson in 
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Figure 7. Magnetization m (cf Eqn. (13)). The error increases significantly as the 
isotropic Heisenberg model is approached. The restoration of the full SU (2) symmetry 
at that point leads to a significant decrease of the magnetization, which is not correctly 
captured by iPEPS for the M we considered. Finite size scaling was performed for the 
Monte Carlo calculations. 



the system. The nature of the transition from the superfluid to the checkerboard-solid 
phase, referred to as spin-flop transition, was under debate for some time until Monte 
Carlo simulations clearly established it as first-order [43j . 

The Hamiltonian has a U{1) symmetry which is spontaneously broken in the 
superfluid phase. At = Jxy^ the symmetry is enlarged to SU{2)^ which is also broken 
with long-range order occurring at T = 0. The enlarged symmetry group leads to larger 
quantum fluctuations, thereby reducing the staggered magnetization as the isotropic 
point is approached. As opposed to calculations on finite systems, the infinite system 
will fully break the symmetry, thereby allowing a direct calculation of the magnetizations 
from local observables. 

The results we obtain for the spin-flop transition are shown in Fig. [5| The results 
for M = 2 already match very well with the QMC results. The change with M = 3 is 
very minor; at /i = 1.5, the relative difference between the critical coupling is on the 
order of 0.1%. A very small M is therefore sufficient to characterize this phase transition 
with very good accuracy. 



3.2. XXX and XY model 

A family of gapless systems is obtained for Jxy = l^h = 0^ ^ [0, 1]. In Fig. |6| results 
for the energy and the spontaneous magnetization are shown, where the magnetization 
m is defined as (subscript indices denote the site in the unit cell) 

^ ^(^/c) unit cell (13) 

i=x,y,z 
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Figure 8. (color online) Different dimerization and unit cell patterns. The pattern 
on the left, which we refer to as staggered pattern, is described by a unit cell with two 
independent tensors and used for the dimerized Heisenberg model. For J' = 0, the 
honeycomb lattice is obtained. The second pattern, referred to as columnar pattern, 
is used in the frustrated cases, but requires 4 instead of two independent tensors. The 
thick, solid bonds carry couplings J^^, J^, the dashed lines carry couplings Jxy^ Jz- 



{a^) is vanishes for all < Jxy. The data clearly shows that while energies are captured 
very well for both models, the large quantum fluctuations away from the Neel order, 
in particular for the isotropic case, are not captured by a PEPS of small dimension 
very well, leading to a significant error in the magnetization. Interestingly, the energy 
for the Heisenberg model decreases strictly monotonously with increasing M, while the 
magnetization does not improve significantly from M = 3 to M = 4 (this can also 
be seen in Fig. [9]). While the discrepancy between the improvement in energy and 
magnetization is clearly due to the vanishing gap above the ground state, the reason for 
the plateau in the magnetization cannot be determined conclusively. 

Figure [7] shows how the decrease of the magnetization as Jz is tuned from to Jxy 
is captured by the PEPS; in particular the rapid decrease as the SU{2) symmetry is 
restored is not seen in the PEPS calculation. 



3.3. Dimerized Heisenberg model 

In order to test the accuracy of iPEPS for a quantum phase transition, we consider a 
dimerized Heisenberg model given by the Hamiltonian 

// = 5]j<^'^V,a„ (15) 

where we choose the couplings J^^'-^^ inhomogenously according to the pattern shown 
on the left in Fig. [S] For J' = 1, J = 0, the state is made up of singlets on the 
strong-coupling bonds. In the limit J = J^, the isotropic Heisenberg model is recovered. 
Between those limits, a second-order phase transition occurs [44j. The behaviour of the 
order parameter, the staggered magnetization, is shown in Fig. |9j While the results do 
indicate a second-order phase transition, the critical point observed for small M deviates 
significantly from the Monte Carlo result. 

This is a clear indication that the phase on both sides of the transition is a highly 
entangled one: even in the limit of pure dimers, M = 2 is required to describe the ground 
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Figure 9. (color online) Staggered magnetization of the dimerized Heisenberg model 
as a function of the coupling ratio. J/J' = l corresponds to the isotropic Heisenberg 
model, while J/ J' = corresponds to isolated dimers on every second horizontal bond. 



state; in the limit of the Heisenberg models, quantum fluctuations around the Neel state 
are very strong. It is interesting to note that the results for M = 3 interpolates between 
those for M = 2 towards the dimerized phase and M = 4 towards the antiferromagnet. 

4. Frustrated Heisenberg Models 

To study the applicability of the iPEPS method to frustrated quantum models, we 
have studied a phase transition in the family of models obtained by choosing the signs 
of the different couplings as shown in Fig. |8] in such a way that the product around 
each plaquette is negative, i.e. JJ^ < 0. The class of models allows a large degree 
of freedom in the choice of parameters. In particular, we can choose full Heisenberg 
coupling Jxy'^^ = Jz'^^ or restrict the coupling to the XY plane, \Jxy\ > 0, = 0. 

The case of Heisenberg couplings on the staggered pattern has been studied by 
Kriiger et al [45] using the Coupled-Cluster Method and Exact Diagonalization. The 
authors find a phase transition from the antiferromagnetic state in the non-frustrated 
limit to a hexatic phase. Classically, this occurs at a coupling ratio J^J = 1/3; in 
the quantum model, the antiferromagnetic phase is stable up to J^J ^ 1.35. While 
the phase transition is second-order in the classical case, it is conjectured to be of 
first order in the quantum case. The hexatic phase for this pattern however has a 
large magnetic unit cell which renders it unsuitable for simulation with iPEPS. In the 
following discussion, we will therefore focus on the case of XY coupling on the columnar 
configuration with J = —1. 

Classical analogues of spin-| systems can be obtained by replacing the classical spins 
with Ising (Z2), XY (0(2)) or Heisenberg (0(3)) spins. Equivalently, one can consider 
mean-field solutions of the quantum system since for spin-| degrees of freedom, the 
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Figure 10. (color online) Energies for the different bonds and comparison of energies 
for the classical and quantum columnar frustrated XY model. 



mean-field approximation to a Hamiltonian of tlie form 

H = J2 J^'^^^j + 4''^^^^] + Jf'^^t^h (16) 

{hj) 

can be mapped exactly to a classical Hamiltonian for Heisenberg spins. We can therefore 
choose to either solve a classical system or perform a mean-field simulation, e.g. by 
minimizing the energy of an iPEPS state with bond dimension M = 1. In this case, the 
renormalization procedure described above reduces to the multiplication with a scalar 
which cancels for all physical observables. 

The classical system has been studied by Villain [46j and solved for the case of 
I J' I = I J| with Ising and XY spins. For XY spins on the staggered pattern, a two- fold 
degenerate hexatic state is found, which exhibits a 2 x 2 site unit cell. Using numerical 
mean-field calculations on a 2 x 2 unit cell, we locate a second-order phase transition 
from a ferromagnetically ordered state to the hexatic state at J^J =1/3. Our result 
for jy J = 1 complies with the result by Villain. 

To simulate the quantum system beyond mean-field, we use an iPEPS ansatz with 
M = 2 and a 2 x 2 unit cell. The classical ground state is used as initial state 
for the iPEPS simulation, since this improves stability and convergence speed of the 
algorithm. We measure the energy on each bond and the relative angles between 
the spins in the XY plane. For the ferromagnetic configuration, these are 0; for 
Z{AC) = Z{AB) = 7r^Z{AD) = 0, the state is antiferromagnetic. In between these 
limits, a hexatic state is found. Therefore, the angles can be interpreted as order 
parameters for the phase transition. 

In Fig. 10, the energy for each bond and a comparison of the total energy per site to 
the classical result is shown. For the non-frustrated case, = 0, we can also compare to 
Quantum Monte Carlo results. The energy found with iPEPS for M = 2 is £"0 = —1.676, 
while the Monte Carlo results for small lattices extrapolate to £"0 = —1.683(1). The 
relative error of AE = 0.0042 is comparable to that obtained for the XY model on the 
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Figure 11. (color online) Relative angles between the magnetizations in the XY plane 
for the quantum and classical columnar frustrated XY model, which is equivalent to 
a mean-field solution of the quantum model. Notice that at maximal frustration, the 
angles for the classical and quantum case coincide. The transition is shifted from 
J' /J = 1/3 to a larger value of J' /J ^ 0.45. 




Figure 12. (color online) Ground state of the columnar-pattern frustrated model with 
Jxy = — 1 and J'^y = 1. Four magnetic unit cells are shown for illustration purposes. 
The arrows indicate the magnetization in the XY plane (up to global rotation). The 
angles of the quantum ground state are identical to those of the classical ground state. 



square lattice. 

We find that the magnetization always remains in the XY plane and is reduced 
by quantum fiuct nations to a value on the order of m 0.45; at all points except 
the maximally frustrated point = J, the magnetization is different on the AD and 



BC sublattices. The angles as a function of the coupling ratio are shown in Fig. 11 



They indicate a phase transition analogous to the classical one, which however is shifted 
towards a higher critical coupling of around / J = 0.45. This agrees with previous 
results indicating that quantum fiuctuations stabilize the collinear order in this class of 
models. 



The state at the point of maximal frustration, \J^/ J\ = 1, is shown in Fig. 12, The 
state corresponds directly to Villain's solution up to a reduction of the magnetization 
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on each site to m ^ 0.451. We have cross-checked this result by starting the simulation 
with different initial states. 

In order to investigate the order of the phase transition, we can study the stability of 
a state upon quenches of the parameters that take the system across the phase transition 
as discussed in Sect. |2.2[ In our simulations, we do not observe such stability of the 
state, indicating that the transition remains continuous in the presence of quantum 
fluctuations. 

5. Conclusion 

We have confirmed that the iPEPS method can be applied with excellent accuracy to 
first-order phase transitions using the example of the spin-fiop transition in the hard- 
core boson model. The stability of a homogeneous system across the phase transition 
allows us to clearly distinguish the nature of the phase transition and determine the 
location of the transition very accurately. 

For the isotropic Heisenberg model and a dimerized Heisenberg model, the accuracy 
compared to finite-size extrapolations of unbiased Monte Carlo simulations decreases 
significantly. This can be understood as a signature of large quantum fiuctuations away 
from the product state, which cannot be captured by a PEPS of small bond dimension. 

For a simple frustrated model on the square lattice, we have shown that the method 
converges to results that agree with expectations based on mean-field simulations and 
previous results obtained using the Coupled-Cluster Method. Whereas the applicability 
of path integral Monte Carlo methods to such systems is limited by the so-called sign 
problem^ the accuracy of iPEPS is only limited by the entanglement of the ground state 
and frustration does not introduce additional difliculties. We therefore expect that the 
method will become a valuable tool for such systems. 

During completion of this work, we became aware of another study of the spin-flop 
transition using an algorithm based on tensor-product states by P. Chen et al [25j . 
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